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Abstract.  The  existing  stability  theory  for  finite  difference  models  of  hyper¬ 
bolic  initial  boundary  value  problems,  due  to  Gustafsson,  Kreiss,  and  Sundstrom, 
is  difficult  to  understand  in  its  original  algebraic  formulation.  Here  we  show 
that  the  GKS  stability  criterion  has  a  physical  interpretation  in  terms  of  group 
velocity:  if  the  finite  difference  model  together  with  its  boundary  conditions  can 
support  a  set  of  waves  at  the  boundary  with  group  velocities  pointing  into  the 
field,  then  it  is  unstable.  A  simple  argument  explains  why  such  a  set  of  waves  is 
unstable,  and  yields  a  new  theorem  on  what  kind  of  unstable  growth  to  expect. 
We  give  examples  in  one  and  two  space  dimensions. 
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0.  INTRODUCTION 


When  time-dependent  partial  dilTerential  equations  arc  solved  numerically  by 
finite  difference  methods,  more  boundary  conditions  are  usually  required  than  the 
physics  of  the  problem  provides.  This  necessitates  a  selection  of  additional  numeri¬ 
cal  boundary  conditions,  and  in  this  choice  an  overriding  consideration  is  numerical 
stability.  For  the  case  of  hyperbolic  equations,  the  stability  question  is  solved  in  prin¬ 
ciple  by  the  theory  of  Gustafsson,  Kreiss,  and  Sundstrom  [9] — henceforth  “GKS”. 
However,  application  of  this  theory  has  been  hindered  by  its  complexity  and  abstract¬ 
ness.  The  purpose  of  this  paper  is  to  point  out  that  the  main  result  of  the  GKS  theory 
has  a  simple  physical  interpretation  in  terms  of  group  velocity.  This  interpretation 
does  not  provide  an  alternative  to  the  algebraic  stability  test  of  the  GKS  theory,  which 
may  in  practice  be  difficult  to  carry  out,  but  it  makes  the  meaning  of  the  algebra  clear. 
It  also  leads  to  a  proof  that  except  for  certain  borderline  cases,  GKS-unstable  models 
are  unstable  not  only  in  the  complicated  norm  of  [9],  but  in  the  ^2  norm  also. 

Group  velocity  is  a  concept  associated  with  energy  propagation  under  dispersive 
equations,  not  hyperbolic  ones.  Its  significance  to  numerical  stability  results  from 
the  fact  that  finite  difference  models,  even  of  nondispersive  equations,  are  necessarily 
dispe.  Ive.  A  general  discussion  of  group  velocity  effects  in  Tinite  difference  models 
may  be  found  in  [13].  For  this  one  should  also  see  the  work  of  R.  Vichnevetsky, 
including  [15]  and  the  papers  referenced  there.  The  particular  topic  of  the  present 
paper  is  treated  in  much  greater  detail  in  the  author’s  PhD  dissertation,  [14]. 

In  the  first  section  here,  we  look  at  a  simple  example  to  illustrate  the  group 
velocity  idea.  Section  2  extends  this  to  the  general  setting  of  GKS.  Section  3  applies 
these  results  to  define  certain  general  classes  of  unstable  difference  models,  in  several 
space  dimensions  as  well  as  one,  which  generalize  examples  that  have  appeared  in  the 
literature. 


1.  AN  EXAMPLE 

Consider  the  model  problem 

Ut=ti*  (l.l) 

for  X  G  (—00,00),  t  G  [0, 00),  with  initial  conditions 

w(x,0)  =  /(x).  (1.2) 

Such  an  initial-value  problem  on  a  domain  without  boundary  is  called  a  Cauchy 
problem.  The  exact  solution  is  u(x^t)  s  f{x  +  ^),  a  leftward  translation  at  speed  1. 


In  a  Fourier  or  normal  mode  analysis  of  (1.1),  one  looks  Tor  wave  solutions 


u(x,  =s  e*^*^*^^*^ 


(1.3) 


where  w  is  the  (temporal)  frequency  and  (  is  the  wave  number.  Eq.  (1.1)  implies  that 
w  and  (  are  related  by 

^  =  (1.4) 

which  is  called  the  dispersion  relation  of  (1.1). 

Let  us  set  up  a  uniform  grid  with  space  step  h  and  time  step  it,  and  approximate 
u(x,  t)  by  a  grid  function  v”: 

vy  »  u{jh^nk)y  (— oo  <  j  <  oo,  n  >  0), 

One  of  the  simplest  difference  models  for  (1.1)  is  the  leap  frog  (LF)  formula 

(1.5) 


where  X  is  the  mesh  ratio  or  Courant  number  k/h.  For  X  <  1  this  scheme  is  stable;  we 
say  that  LF  is  Cauchy  stable.  By  plugging  (1.3)  into  (1.5),  one  obtains  the  dispersion 
relation  for  LF, 

=  e-*"*'  + 


that  is, 


sin  wk  =  —X  sin  ^h. 


(1.6) 


This  relation  approximates  (1.4)  only  for  small  wk  and  ^h — that  is,  only  for  waves 
that  are  well  resolved  on  the  grid. 

Since  w  is  no  longer  a  linear  function  of  f ,  the  LF  model  is  said  to  be  dispersive. 
According  to  a  theory  going  back  to  Lord  Rayleigh  [3,4,11,16],  energy  associated  with 
wave  number  (  will  propagate  at  a  group  velocity*  given  by 


C=: 


duj 


(1.7) 


For  LF  one  obtains,  by  differentiating  (1.6)  implicitly, 

(1.8) 

It  is  apparent  that  energy  associated  with  different  wave  numbers  or  frequencies  will 
travel  at  different  group  velocities,  so  that  an  initial  signal  that  is  not  monochromatic 
will  change  form  as  it  propagates.  For  a  well  resolved  wave  one  has  «  0  and 

*The  group  velocity  it  not  the  tame  at  the  phate  velocity,  C  =  Rcadcrt  unramtitar  with  this 

diatinction  are  encouraged  to  read  the  ditcuttiont  in  [4],  [11],  or  [16]. 


COS  ^h 
coswib* 
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(jjk  0,  hence  wfc  «  — X^/i  by  (1.6),  and  (1.8)  becomes 

c  =  -1  +  +  o{uh)*).  (1.9) 

Thus  typical  signals  travel  leftward  at  a  speed  less  than  the  ideal  speed  1,  with 
higher  wave  numbers  lagging  more  than  lower  ones.  This  dispersion  of  wave  numbers 
gives  rise  to  the  spurious  oscillations  near  discontinuities  that  arc  familiar  in  finite 
difference  computations.  Alternately,  if  one  sets  up  a  wave  packet  consisting  of  energy 
at  essentially  constant  w  and  the  packet  will  be  seen  to  move  leftward  without 
changing  at  the  velocity  (1.8)j=»(1.9).  For  illustrations  see  [13,15]. 

The  key  to  our  analysis  is  that  (1.8)  is  valid  not  only  for  well  resolved  waves,  but 
for  all  waves  supportable  on  the  grid,  which  means  a  range  of  {— 7r,7r]  in  both  and 
(jjk.  In  fact  for  many  waves,  C  has  the  wrong  sign,  so  that  energy  travels  in  the  wrong 
direction.  In  particuar,  (1.8)  gives  the  following  group  velocities  for  four  extreme 
situations — the  constant  function  and  three  parasitic  waves  that  are  sawtoothed  in  x 
and/or  t: 

Uh,wk) 

(а)  (0,0)  -1 

(б)  (tt.O)  1 

(c)  (0,x)  1 

(d)  (-ir,ir)  -1 


Line  (b)  implies,  for  example,  that  if  initial  data  are  supplied  to  LF  of  the  form 


-1)^'  [jh  <  c) 

0  (jh  >  e) 


(1.10) 


for  some  c  ^  A  >  0,  then  the  result  as  t  increases  will  be  a  steady  rightward  motion 
of  the  wave  front  from  x  =  c  at  speed  1.  See  Fig.  6  of  [13j. 

Now  let  us  turn  to  an  initial  boundary  value  problem.  Let  (1.1)  and  (1.2)  be 
given  on  the  quarter-plane  >  0;  no  boundary  data  at  x  =  0  are  needed  to  make 
this  problem  well  posed.  To  obtain  an  approximate  solution  on  the  gnd  j,n  >  0, 
wc  can  specify  initial  values  and  Vy  for  j  >  0,  and  apply  LF  for  n  >  2  at  points 
j  >  1.  An  additional  boundary  formula  is  then  needed  for  Uq,  n  >  2. 

Suppose  we  (foolishly)  pick  the  bound.ary  formula 


=  +  (n>l) 


(1.11) 


and  proceed  to  step  forward  in  time.  Now  imagine  that  at  some  time  step  a  pertur¬ 
bation  (e.g.  rounding  error)  happens  to  be  introduced  that  has  the  form  (1.10).  It  is 
easy  to  sec  that  such  a  wave  satisfies  (1.11),  and  therefore  this  mode  will  behave  just 
as  if  the  domain  were  still  (—00,  00):  the  wave  front  will  begin  to  propagate  rightwards 
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n  *  1 


t  «  .01?.5 


n  »  4 
t  *  .05 


n  *  49 
t  =  .5 


FIG  1.  Appearance  of  an  unstable  rightgoing  parasitic  wave 

in  the  Leap  Frog  model  (1.5),  (1.11)  with  initial  conditions  (1.12). 


into  r  >  0  at  speed  1.  The  initial  perturbation,  with  sum-of-squares  energy  on  the 
order  of  €,  will  give  rise  to  a  growing  solution  with  energy  on  the  order  of  e  + 
Since  e  might  be  arbitrarily  small,  this  amounts  to  an  amplification  of  the  initial 
p<^rturbation  by  an  unbounded  factor.  The  difference  scheme  is  unstable,  because 
there  exists  a  rightgoing  wave  that  satisifies  both  the  interior  formula  (1.5)  and  the 
boundary  condition  (l.ll)* 

One  can  verify  experimentally  that  the  scheme  (1.5),  (1.11)  is  susceptible  to  an 
unstable  rightgoing  mode  of  type  (b).  Fig.  1  shows  a  computation  on  a  grid  with 
h  =  1/‘10,  X  =  1/2.  For  initial  data  we  took  Oy  =  vj  =  0  for  all  j  except  for  the 
“random'*  nonzero  initial  values 

=  1.  ^0  =  (1.12) 

Figs,  la-c  show  the  resulting  solution  at  steps  n  =  1,4,40,  i.e.  t  =  .0125,  .05,  .5. 
Obviously  the  expected  mode  has  been  excited,  and  apparently  no  others.  In  a  realistic 
computation,  rounding  errors  would  cause  a  similar  radiation  of  energy  in  this  mode 
from  the  boundary.  From  (1.8)  one  can  see  that  there  are  many  other  rightgoing 
modes  for  LF — in  fact,  any  wave  with  ^h  <  tt/2  and  wk  >  7r/2  or  ^h  >  7r/2  and 
ufk  <  7r/2.  Mode  (c)  is  the  simplest  example.  None  of  these  lead  to  instabilities, 
however,  because  none  of  them  satisfy  (1.11). 

The  main  GKS  theorem  asserts,  roughly,  that  an  initial  boundary  value  problem 
model  is  stable  if  and  only  if 

(i)  the  interior  difference  formula  is  Cauchy  stable; 

(ii)  the  model  (including  boundary  conditions)  admits  no  eigcnsoluUons  that  grow 
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from  each  time  step  to  the  next  by  a  constant  factor  z  with  \z\  >  1; 

(iii)  the  model  (including  boundary  conditions)  admits  no  wavelike  solutions  with 
group  velocity  C*  >  0. 

In  practice,  (i)  is  easy  to  verify,  and  there  are  rarely  any  growing  modes  of  type  (ii). 
The  critical  condition  is  therefore  (iii).  We  will  state  the  GKS  theory  more  precisely 
in  the  next  section* 

To  obtain  stability  in  the  present  example,  we  might  replace  (1.11)  by  the  con¬ 
dition 

Then  it  is  not  hard  to  verify  that  (ii)  is  satisfied,  and  the  only  wavelike  mode  admitted 
by  the  difference  model  is  (^h,  wfc)  =  (0,0),  which  has  C  =  <  0,  hence  cannot 

cause  radiation  from  the  boundary.  Thus  (iii)  is  satisfied  and  the  model  is  stable. 

It  may  seem  that  the  kind  of  energy  growth  in  Fig.  1  is  too  weak  to  be  dangerous, 
and  should  not  be  considered  unstable.  Though  this  may  be  true  for  the  present  model 
problem,  it  is  usually  not  true  in  realistic  computations,  for  various  reasons.  One  is 
that  when  a  rightgoing  signal  of  constant  amplitude  can  occur,  then  sometimes  similar 
signals  with  amplitudes  growing  linearly  with  n  or  faster  are  also  possible.  Such  noise 
might  soon  become  significant  even  though  it  began  at  the  level  of  rounding  error. 
Second,  if  the  semiinfinite  region  x  >  0  is  replaced  by  a  bounded  strip  such  as  (0, 1), 
then  a  wavelike  instability  of  type  (iii)  may  be  converted  by  repeated  reflections  back 
and  forth  to  an  exponential  growth  of  type  (ii),  which  is  unambiguously  unstable. 
Third,  the  rate  of  growth  may  also  be  increased  if  variable  coefficients  or  lower  order 
(undifferentiated)  terms  are  present.  The  GKS  theory  shows  that  by  ruling  out  the 
relatively  weak  instability  of  Fig.  1,  one  can  be  certain  that  these  more  serious  effects 
arc  also  excluded. 

In  summary,  GKS  instability  amounts  to  the  spontaneous  radiation  of  energy 
from  the  boundary  into  the  interior.  In  a  realistic  physical  application,  this  is  likely 
to  cause  trouble. 


2.  GENERAL  STATEMENT  OF  THE  GKS  THEORY 


Consider  now  a  first-order  hyperbolic  system 

—u{x,t)=z  A—uix,t),  (2.1) 

where  u  is  an  A^-vector  and  >l  is  a  constant  nonsingular  Hermitian  matrix  of  dimension 
N  X  N.  (All  of  what  follows  extends  to  equations  that  include  a  zeroth-order  term 
Bu[x,i)  and  a  forcing  function  F{x,^).)  Let  (2.1)  be  modeled  in  x  >  0  by  a  fixed 
^-h2-leve]  Cauchy  stable  diiference  formula,  explicit  or  implicit,  whose  stencil  extends 
t  points  to  the  left  and  r  points  to  the  right  of  center.  For  example,  LF  has  8=1 
and  /  =  r  =  1.  We  can  write  the  difference  model  formally  as 

a 

=  53  Qcv’'-^  (2.2) 

<F=0 

where  vj  «  w(jh,  nk)  is  an  AT- vector  for  each  j  and  n,  and  each  is  a  difference 
operator  with  matrix  coefficients  of  size  IV  X  N.  For  simplicity  we  will  assume  that 
each  Qa  is  constant,  independent  of  h.  If  (2.2)  is  applied  for  j  >  t,  then  boundary 
formulas  are  required  to  determine  values  for  j  =  0,  1: 

0  <  J  <  <  - 1-  (2.3) 

These  formulas  comprise  both  physical  boundary  conditions,  such  as  couplings  be¬ 
tween  outflowing  and  inflowing  components  of  v,  and  additional  purely  numerical 
boundary  conditions.  Unfortunately,  it  would  take  many  pages  to  state  precisely  the 
form  of  the  difference  model  and  the  assumptions  it  must  satisfy,  so  for  details  the 
reader  is  referred  to  §1  and  §5  of  [9]  and  to  §i  of  (6],  where  the  presentation  is  clearer. 

The  GKS  idea  is  to  perform  a  normal  mode  analysis  of  what  solutions  the  model 
(2.2),  (2.3)  can  support.  To  find  the  normal  modes,  let  us  begin  by  ignoring  the 
boundary  conditions  (2.3).  In  the  last  section,  we  looked  for  wavelike  solutions  (1.3) 
with  frequency  w  and  wave  number  Let  us  now  write 

z  =  =  (2.4) 

The  vector  analog  of  (L3)  then  takes  the  form 

Vy  =  (2.5) 

where  is  a  constant  nonzero  7V-vector.  Now  suppose  that  instead  of  taking  \z\  =  1, 
we  let  z  be  any  complex  number  with  |x|  >  1.  Then  there  are  still  modes  of  the  form 
(2.5),  but  now  k  will  become  complex  also,  with  1.  (The  circumstance  |2r|  >  1, 

\k\  =:  1  would  violate  the  von  Neumann  condition,  and  we  have  assumed  that  (2.2) 


FIG  2.  Sketch  of  signals  with  (zj  >1 


(a)  Rightgoing,  jic|  <1 


(b)  Leftgoing,  1k1  >  1  . 


is  Cauchy  stable.)  Figure  2  suggests  the  two  possibilities  |/c|  <  I  and  |k|  >  1.  From 
one  time  step  to  the  next,  each  mode  in  Figure  2  increases  in  amplitude  by  a  ratio 
|^(.  It  is  obvious,  however,  that  we  may  view  this  equivalently  as  a  lateral  motion, 
rightward  in  case  (a)  and  leftward  in  case  (b),  combined  with  a  change  of  phase.  We 
will  call  these  Tigktgoing  and  leftgoing  modes,  respectively. 

Modes  of  the  form  (2.5)  do  not  always  span  the  set  of  solutions  with  time 
dependence  z”.  If  p  >  1  /c*s  coalesce  for  some  z,  then  a  defective  situation  results  in 
which  a  p-parameter  collection  of  modes  are  possible, 

v;  =  (0  <  d  <  P  -  1).  (2.6) 

(The  defective  situation  p  >  1  rarely  occurs  in  practice.)  A  fundamental  lemma  of 
(9]  now  states 

Proposition  (f9],  Lemma  5.2).  For  any  z  with  |z|  >  1,  (2.2)  admits  a  total  of 
exactly  Nl  linearly  independent  solutions  of  type  (2.6)  with  |/c|  <  1,  and  exactly  Nr 
such  solutions  with  |/c|  >  1.  | 

This  result  is  proved  by  reducing  the  di (Terence  model  to  a  recurrence  relation  in  j. 
From  it  we  get  the  general  solution  to  (2.2)  that  varies  with  z": 

Proposition.  For  any  z  with  |z|  >  1,  the  general  solution  to  (2.2)  of  the  form 
=  2^<l>j  may  6c  written 

=  z:"  ^  ^  QK.dK^j'^rl’K  +  ^  (2.7) 

|/e(»)|<l  4=0  I*(*)I>1  ^<“0 
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for  arbitrary  constants  | 


We  may  think  of  (2.7)  as  describing  a  linear  combination  of  jVf  rightgoing  and  Nr 
leftgoing  modes. 

Now  U:t  us  reintroduce  the  boundary  conditions.  The  full  model  can  support 
only  any  solutions  (2.7)  that  satisfy  (2.3)  as  well  as  (2.2).  There  are  Nt  boundary 
formulas,  and  (2.7)  depends  on  N[t  +  r)  free  parameters,  so  one  may  expect  that 
there  will  normally  be  an  yVr-dimensional  space  of  solutions  that  the  model  supports. 
In  such  a  solution,  a  combination  of  leftgoing  waves  hits  the  boundary  and  reflects 
as  a  combination  of  rightgoing  waves.  This  kind  of  configuration  is  in  general  not 
unstable,  even  though  the  reflected  energy  may  be  larger  than  the  incident  energy. 
But  occasionally  it  may  happen  that  some  collection  of  rightgoing  signals  satisfies 
(2.2)  and  (2.3)  by  itself,  i.e.  without  any  stimulation  by  leftgoing  signals.  In  this 
case  the  left/right  reflection  coefficient  has  become  infinite.  Then  the  model  has  an 
instability  of  the  Godunov- Ryabenkii  (G-R)  type  [12]: 

Godunov- Ryabenkii  theorem.  Suppose  that  for  some  z  with  (2:|  >  1,  the 
model  (2,2),  (2.3)  admits  a  solution  of  the  form 

P«-i 

z”<j>j  =  z"  ^  ^  (2.8) 

|<c|<l  d=0 


Then  it  is  unstable. 

Proof.  If  Vy  =  is  taken  as  initial  data  for  0  <  cr  <  s,  the  solution  as  n 
increases  will  be  for  all  n.  Since  t  =  nk,  this  means  that  v  will  grow  like 

\zY^^.  With  \z\  >  I  this  growth  is  unbounded  for  any  f  as  fc  — ^  0.  | 

A  similar  solution  made  up  of  leftgoing  waves  would  of  course  also  grow  unboundedly 
in  amplitude,  but  this  would  amount  to  propagation  of  energy  leftward  from  infinity 
rather  than  generation  of  energy  at  the  boundary.  The  initial  wave  would  not  have 
finite  size  in  or  any  other  reasonable  norm,  and  could  not  be  excited  by  rounding 
errors  or  other  perturbations.  The  wave  ^  in  (2.8),  on  the  other  hand,  docs  belong  to 
^2,  and  is  called  an  eigensolution  of  the  difference  model. 

The  set  of  potential  eigensolutions  since  they  are  made  up  of  rightgoing 
signals  only,  is  spanned  by  Nt  free  parameters,  which  exactly  matches  the  number  of 
boundary  conditions  (2.3).  Therefore  the  G-R  condition  has  this  algebraic  form:  does 
there  exist,  for  any  >  1,  a  nontrivial  solution  of  a  system  of  Nf  linear  homogeneous 
equations  in  unknowns  that  depends  on  Hence  the  problem  can  be  cast  as  a 
determinant  condition  involving  an  Nt  X  Nt  matrix  M{z).  (For  examples  see  [5,6,9].) 

Godunov-Ryabenkii  theorem  (determinant  condition).  A  necessary  con¬ 
dition  for  stability  of  the  model  (2.2),  (2.3)  is 


det  M[z)  0 
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(2.9) 


for  all  z  with  \z\  >  1.  | 


A  virtue  of  normal  mode  analysis  is  that  in  principle,  the  determinant  condition  can  be 
verified  mechanically,  although  in  practice  this  may  be  very  difficult  [6],  In  contrast, 
stability  proofs  by  the  energy  method  [10,12]  require  a  measure  of  intuition  and  luck. 

The  limitation  of  the  G-R  condition  is  that  it  is  a  necessary  condition  for  stability, 
but  far  from  sufficient.  To  obtain  a  condition  that  is  sufficient  or  nearly  so,  one 
must  investigate  the  borderline  case  |^|  =  1.  This  is  the  main  achievement  of  the 
GKS  theory.  From  the  GKS  point  of  view,  this  investigation  is  a  matter  of  showing 
algebraically  that  certain  resolvent  estimates  obtained  for  |2r|  >  1  extend  continuously 
to  \z\  =  1.  From  our  point  of  view,  it  is  a  matter  of  observing  that  the  ideas  of 
“leftgoing”  and  “rightgoing”  still  make  sense  as  (zj  — ►  1,  because  in  the  limit  \2\  =  1 
the  lateral  motion  in  Figure  2  becomes  a  group  velocity. 

First,  the  GKS  approach.  For  \z\  =  1,  what  unstable  solutions  along  the  lines 
of  (2.7),  but  also  including  components  with  |/c|  =  1,  can  (2.2),  (2.3)  support?  The 
GKS  theory  gives  the  following  answer  based  on  a  perturbation  teat  First,  one  can 
rule  out  components  with  \k\  >  1,  as  before.  One  then  investigates:  are  there  any 
solutions  to  (2.2)  and  (2.3)  of  the  form  (cf.  (2.8)) 

=  ^"  XI  XZ  {2-10) 

|«|<1  d=Q 

for  some  z  with  |^|  >  1?  If  not,  the  model  is  stable.  If  so,  a  further  test  is  required. 
Given  a  solution  (2.10),  let  z  be  perturbed  slightly  to  a  new  value  z  with  \z\  >  1. 
Each  K  with  |«|  =  1  will  then  move  to  a  nearby  value  k  with  \k\  ^  1.  If  every  k 
for  which  7^  0  in  (2.10)  moves  to  \k\  <  1,  then  (2.10)  is  an  unstable  mode  and 
the  model  Is  unstable.  A  solution  <j>  of  this  sort  with  |^|  =  1  and  |/c|  =  1  for  at  least 
one  K  is  called  a  generalized  eigenaolution  of  the  difference  model.  The  qualification 
“generalized”  has  been  introduced  because  such  a  solution  no  longer  belongs  to  ^2* 

The  GKS  stability  theorem  can  now  be  stated: 

GKS  stability  theorem.  The  model  (2.2),  (2.3)  ia  atable  if  and  only  if  (2.2) 
and  (2.3)  admit  no  eigenaolutiona  or  generalized  eigenaolutiona  with  \z\  ^  I-  t 

A  simpler  but  less  constructive  way  to  express  the  same  theorem  is  by  means  of  the 
determinant  condition.  The  Ni  X  Nt  matrix  M{z),  which  for  \z\  >  1  embodies  the 
condition  that  there  be  an  cigensolution  of  G-R  type,  has  a  continuous  extension  to 
\z\  =  1.  Let  M(z)  for  \2\  >  1  denote  this  extension.  Then  one  has: 

GKS  stability  theorem  (determinant  condition.)  A  neceaaary  and  sufficient 
condition  for  atability  of  the  model  (2.2),  (2.3)  ia 

dciM(z}^0 


for  all  z  with  |-?|  >  1. 
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Proof.  This  result  is  stated  as  Lemma  10.3  and  the  sentence  following  in  [9).  | 


Now  for  the  group  velocity  interpretation.  Assume,  as  is  usually  the  case,  that  k 
depends  analytically  on  z  in  a  neighborhood  of  a  point  z  where  |/c(2)|  =  |2r|  =  1.  By 
(2.4),  we  have  then 

dz  =  ikz  dojf  dk  =  —ihn  d^ 

and  therefore 

_  dz/ikz  ^  1  dz/z 

d^  dkjihK  Xdkfn 

If  C  >  0,  (2.11)  implies  that  dzjz  and  d/c/zc  will  have  equal  and  opposite  signs.  In 
other  words,  k  moves  inside  the  unit  circle  when  z  moves  outside.  Conversely,  if 
C  <  0,  zc  moves  outside  the  circle  when  z  does,  and  the  solution  (2.10)  does  not 
represent  a  generalized  eigensolution.  Thus  the  GKS  perturbation  test  for  generalized 
eigensolutions  corresponds  to  a  test  for  positive  group  velocities. 

The  argument  of  the  example  in  Section  1  explains  why  a  solution  made  of  waves 
with  positive  group  velocities  should  be  unstable.  The  same  idea  leads  naturally  to 
the  following  necessary  condition  for  stability  in  the  general  case: 

Theorem  1#  Suppose  that  the  model  (2.2),  (2.3)  admits  a  solution  of  the  form 
(2.10)  with  \z\  =  1,  where  for  each  k  with  [zcj  =  1,  the  group  speed  defined  by  (2.11) 
exists  and  satisfies  >  0.  Suppose  that  for  at  least  one  such  /c,  say  zcq,  >  0. 
Then  the  model  is  unstable.  | 

In  fact,  we  can  be  more  precise: 

Theorem  2.  Suppose  that  the  model  (2.2),  (2.3)  admits  an  unstable  generalized 
eigensolution  as  described  in  Thm.  1.  Then  there  exists  a  constant  />  >  0  such  that 
j/  6  >  0  and  t  '>  0  are  arbitrary,  then  for  all  sufficiently  small  /i,  k,  there  exists  a  set 
of  intitial  data  0  <  a  <  s,  such  that 

I  for  j  >0f  0  <  cr  <  s, 

v?  =  0  for  jh  >  €f  0  <  cr  <  8, 

and 

>  pt,  (2.12) 

where  ||  )|2  denotes  the  norm  defined  by 


J=0 

Proof.  Assume  d  =  0  for  all  waves  in  (2.10).  Let  initial  data  be  chosen  equal 
to  l/a^iax  times  (2.10)  near  x  =  0,  but  cutting  off  smoothly  to  0  for  jh  >  e.  As 
time  elapses,  each  smooth  wave  front  will  propagate  rightwards  at  its  group  speed 

i 


C^.  This  will  cause  growth  at  least  as  fast  as  (2.12)  if  p  is  taken  as 
For  a  more  rigorous  proof,  including  the  case  d  >  0,  see  [14].  J 

Theorems  1  and  2  are  not  quite  the  same  as  the  GKS  stability  theorem.  The 
latter  asserts  that  even  if  every  wave  in  (2.10)  has  C  =  0,  the  model  is  still  unstable. 
From  the  group  velocity  point  of  view,  it  is  not  clear  why  this  should  be  so,  for  if  a  cut 
off  wave  is  set  up  as  in  the  proof  of  Thm.  2,  but  with  C  =  0  for  each  component,  then 
as  t  increases  the  wave  front  will  approximately  remain  stationary.  The  explanation 
is  that  the  GKS  theorem  does  not  define  stability  in  terms  of  a  simple  I2  norm,  but 
in  a  more  complicated  fashion  under  which  a  wave  that  remains  stationary  at  the 
boundary  turns  out  to  be  unstable  (Defn.  3.3  of  [9]).  The  GKS  definition  has  other 
peculiarities  too,  the  most  troublesome  of  which  is  that  it  defines  stability  not  in 
terms  of  a  norm  at  a  fixed  time  /,  but  in  terms  of  an  integral  of  the  solution  over  all 
t  >  0,  We  will  not  go  into  the  details. 

The  situation  here  is  something  like  that  discussed  at  the  end  of  Section  1:  by 
choosing  a  conservative  definition  of  stability,  one  can  obtain  a  theory  that  is  robust 
with  respect  to  variable  coefficients,  undifferentiated  terms,  and  so  on.  In  fact,  the 
striking  difference  between  the  GKS  theorem  and  our  Thm.  1  (coupled  with  the  G-R 
condition)  is  that  the  former  gives  a  necessary  and  sufficient  condition  for  stability, 
made  possible  by  the  use  of  the  complicated  norm.  The  difference  between  this 
situation  and  that  of  Section  1  is  that  in  this  case  it  is  not  as  clear  that  in  realistic 
computations,  the  conservative  choice  is  generally  necessary.  We  are  now  dealing  with 
a  borderline  subcase  (C7  =  0)  of  a  borderline  case  (wavelike  eigensolutions).  There  do 
not  appear  to  be  any  published  examples  of  instabilities  of  the  0  =  0  type  in  which 
numerical  trouble  actually  becomes  evident. 

In  practice,  stability  for  a  particular  problem  will  be  determined  by  a  host  of  in¬ 
teracting  phenomena  that  depend  on  whether  and  how  smoothly  the  coefficients  vary, 
whether  an  undifferentiated  term  is  present,  whether  homogeneous  or  inhomogeneous 
boundary  data  are  supplied,  whether  there  is  more  than  one  boundary,  and  whether 
the  boundaries  are  characteristic.  Probably  no  simple  theory  can  encompass  all  per¬ 
mutations  of  such  effects  without  some  artificiality.  For  a  further  discussion  of  these 
matters,  see  [14]. 
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3.  APPLICATIONS 


We  will  now  give  five  examples  of  unstable  boundary  conditions,  or  their  equiv¬ 
alents.  In  each  case  the  unstable  mode  consists  of  nothing  more  than  a  combination 
of  constant  and  sawtoothed  signals,  and  this  is  typical  for  instabilities  encountered  in 
practice.  In  the  course  of  the  examples  we  will  extend  the  theory  to  problems  with 
interfaces  and  to  multidimensional  problems,  where  a  vector  group  velocity  becomes 
needed. 

Because  of  the  repeated  occurence  of  sawtoothed  waves,  it  is  convenient  to  devise 
a  name  for  schemes  under  which  they  propagate  in  the  wrong  direction: 

Defn.  Let  Q  be  a  scalar  difference  formula.  Suppose  that  whenever  Q  admits 
a  solution  with  \z\  =  I  and  group  speed  C  €  then  it  also  admits  the 

solution  v"  =  and  this  wave  has  group  speed  C"  G  K  satisfying  CC^  <  0. 

Then  Q  is  x-reversing.  Likewise  if  the  existence  of  a  solution  v"  =  with  |/c|  =  1 
and  group  speed  C  implies  the  existence  of  a  solution  v”  =  (— with  CC*  <  0, 
then  Q  is  t-reversing. 

Any  scheme  based  on  the  standard  second- order  centered  difference  in  x  or  ^  is 
reversing  for  that  variable,  with  CC*  ^  —1.  Thus  LF,  CN  (Crank- Nicolson),  and 
BE  (Backwards  Euler)  are  x- reversing,  and  LF  and  LF4  (fourth-order  leap  frog) 
are  /-reversing,  as  is  any  modification  of  LF  with  (spatial)  dissipation  added.  More 
generally,  the  general  (2/-*- 1)- point  difference  approximation  to  d/dx  or  d/dt  of  order 
2t  is  also  reversing  [14],  with  CC'  <  —  I  for  ^  >  2.  A  dissipative  scheme  cannot  be 
X- reversing,  and  a  scheme  that  dissipates  oscillations  in  t  (for  which  there  is  no  name) 
cannot  be  /-reversing  |14). 

Application  1:  space  extrapolation  with  /-reversing  formulas 

Let  (1.1)  be  modeled  by  a  difference  formula  Q  for  j  ^  t  coupled  with  qjth-ordtr 
apace  extrapolation  boundary  conditions 

S  :  [{E  -  =0  (0  <  i  <  /  -  1) 

for  the  boundary  points,  where  E  is  the  shift  operator  defined  by  [Ev]^  =  and 
qj  >  I  for  each  The  result  appears  in  [9],  and  in  various  other  papers,  that  S  is 
unstable  if  £  =  1  and  the  interior  scheme  is  LF.  Here  is  a  generalization: 

Theorem  3.  Any  eonaistent  t-reversing  difference  formula  Q  for  (1^1),  such  as 
LF  or  LF4  u)ith  or  without  disaipation,  is  unstable  m  combination  with  the  boundary 
condition  S. 

Proof  The  sawtoothed  wave  t;”  =  (“!)’'  satisfies  S  for  any  set  {qj},  and  if  Q 
is  /-reversing,  it  also  satisfies  Q  and  has  C  >  0,  since  by  consistency  v”  =  1  must 
satisfy  Q  with  C  =  —  1  <  0.  By  Thm.  1,  the  model  is  therefore  unstable.  | 


The  instability  with  S  of  an  LF  scheme  with  spatial  dissipation  added  is  pointed  out 
by  Goldberg  and  Tadmor  in  (8). 


Application  2:  ^one-sided  leap  frog”  with  ^^reversing  formulas 

Similarly,  it  has  been  noted  in  various  papers  that  if  (1.1)  is  modeled  by  LF  for 
y  >  1  with  the  boundary  condition 


t,5+‘  =  wr  ‘  +  2X(t;?  -  rj), 


then  the  result  is  GKS-unstable.  As  a  generalization,  consider  any  set  of  boundary 
conditions 


+  ikDjv'^  0  <  i  <  /  -  1,  (3.1) 

where  each  Dj  is  a  one-sided  spatial  difference  operator  involving  at  most  j  points  to 
the  left  of  center,  consistent  with  d/dx.  We  obtain  just  as  above 

Theorem  4.  Any  consistent  t-reversing  difference  formula  Q  for  (1.1)  is  unstable 
in  combination  with  the  boundary  condition  (S.l).  | 


Application  3:  sign-changing  coefficients;  nonlinear  instability 
Consider  the  problem 


f  a-ti*  (x  <  0) 

\  O+Tix  (x  >  0), 


(3.2) 


where  o«  and  are  constants.  To  model  this  by  finite  differences,  we  might  set 
up  a  grid  ((j  -I-  |)fi,  nfc)  for  — oo  <  j  <  oo,  n  >  0.  Suppose  we  apply  consistent 

difference  formulas  Q_  and  for  x  <  —h/2  and  x  >  /i/2,  respecMvely,  taking 

no  special  measures  to  improve  accuracy  at  the  interface.  The  stability  question  for 

such  an  interface  is  essentially  the  same  as  for  an  initial  boundary  value  problem, 

and  the  GKS  theory  has  been  applied  to  such  problems  by  Kreiss,  Ciment,  Oliger, 
and  others  [5].  Formally,  a  scalar  model  including  an  interface  can  be  “folded”  into 
an  initial  boundary-value  problem  for  a  system  of  two  variables,  and  then  the  GKS 
theory  is  directly  applicable.  What  is  really  going  on  in  such  a  process  is  a  search  for 
eigensolutions  or  generalized  eigensolutions  consisting  of  waves  that  are  outgoing  from 
the  point  of  view  of  the  interface.  That  is,  a  difference  model  involving  a  scheme- 
or  mesh-change  interface  is  GKS-stable  if  and  only  if  there  is  no  eigensolution  of  the 
kind  suggested  in  Fig,  S: 


^  FIG  3.  In  a  problem  containing  an 

4 - -  interface,  a  solution  consisting  of 

^  a  set  of  outgoing  waves  on  each  side 

^  will  be  unstable. 


If  sgna_  =  sgna^.,  then  most  models  for  (3.2)  are  stable,  but  stability  vanishes 
if  sgna-  ^  sgna+: 

Theorem  5.  Let  (3,2)  be  modeled  by  consistent  formulas  and  as  indi¬ 
cated  a6ot7C.  //a—  >  0  >  the  model  is  unstable.  If  a  ^  <  0  <  and  Q—  and 
are  both  x -reversing  or  both  t- reversing,  the  model  is  also  unstable. 

Proof  In  the  first  case,  the  constant  function  =  1  is  an  outgoing  wave  that 
satisfies  all  the  difference  formulas,  so  the  model  is  unstable  by  Thm.  1.  In  the  second 
case,  the  same  goes  for  a  space  or  time  sawtooth  (—1)^  or  (—1)".  | 

This  elementary  example  is  related  to  certain  known  examples  of  nonlinear 
instability.  If  the  Burgers  equation 


Ut  =  utt* 


is  modeled  by  the  LF  scheme 

vy  “  l^y+i  “ 

then  exponentially  growing  instabilities  arise  that  are  characterized  by  oscillations  of 
the  form  [7,10] 

^  0,  <  0,  ^7+3  ^  0. 

Though  it  is  easy  to  enough  to  examine  this  problem  directly,  it  also  has  a  rough 
interpretation  along  GKS  lines.  LF  is  an  x-reversing  formula,  and  the  instability 
observed  looks  approximately  like  the  outgoing  sawtooth  of  Thm.  5  from  the  point  of 
view  of  the  sign-change  interface  at  Xj>3/?.  The  linear  growth  of  this  outgoing  wave 
would  be  converted  into  exponential  by  reflection  at  points  Xj  and  xy+3  even  if  the 
coefficients  Vj  did  not  change  from  one  time  step  to  the  next;  the  fact  that  they  do 
makes  the  growth  even  more  rapid. 
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Application  4:  mesh  refinement 

In  problems  where  the  solution  is  smoother  in  some  regions  than  in  others,  it 
may  be  useful  to  combine  grids  of  various  sizes  [5].  Where  two  such  grids  meet,  some 
kind  of  interface  condition  will  be  required.  Here  is  a  mesh  refinement  scheme  for 
which  the  stability  result  takes  a  particularly  interesting  form  (Joseph  Oliger,  private 
communication).  Suppose  that  for  x  >  0,  the  grid  Xj  =  hj  {j  >  0)  is  set  up,  while  for 
X  <  0,  this  is  coarsened  by  an  integral  factor  of  m  >  2,  so  that  the  grid  is  xy  =  mhj 
for  j  <  0.  If  (1.1)  is  modeled  by  LF  for  j  <  and  j  >  If  a  formula  is  still  needed 
to  determine  tiq.  The  obvious  choice  is  to  apply  the  coarse  grid  LF  formula  at  j  =  0, 
which  Is  possible  because  the  value  required  at  x  =  m/i  is  available  from  the  fine  grid: 

v5+‘=t,5-i+mX(C-v!l,).  (3.3) 


Now  suppose  a  wave  is  considered  of  the  form 


(i>0), 
1  1  (i<o). 


(3.4) 


On  X  <  0,  this  wave  is  constant  and  has  C  =  —1.  On  x  >  0,  it  is  sawtoothed  and  has 
C  >  0  since  LF  is  x- reversing.  Thus  (3.4)  is  outgoing  on  both  sides  of  the  interface. 
Moreover  if  m  is  even,  it  obviously  satisfies  (3.3),  so  we  have  instability. 

This  conclusion  can  be  generalized  as  follows: 

Theorem  6.  Let  (1>1)  be  modeled  by  a  consistent  x-reversing  S-point  formula 
on  Xj  =  jh  for  J  >  1  coupled  with  any  consistent  formula  on  xy  =  jmh  for  j  <  0, 
with  right-hand  values  for  the  latter  near  the  interface  taken  where  needed  from  points 
imh  with  i  >  I.  Then  if  m  is  even,  the  model  is  unstable*  | 

For  LF  or  CN  the  sawtooth  (3.4)  turns  out  to  be  the  only  instability  that  arises,  so 
this  kind  of  mesh  refinement  is  stable  if  m  is  odd. 


Application  5:  two^^dimensional  problems 

Abarbanel  and  Gottlieb  [1]  and  Abarbanel  and  Murman  [2]  have  studied  the 
stability  of  various  dilTcrence  schemes  for  the  following  problem  in  two  space  dimen¬ 
sions: 

Ut  u,  +  x,t  >  0,  1/6  (—00,00).  (3.5) 

The  solutions  to  this  equation  consist  of  functions 

u(i,  y,t)  =  u(x  +  t,y  +  t,0). 

That  is,  information  propagates  with  a  vector  velocity  (—1,-1).  Since  the  (low  is 
outward  across  the  boundary  z  =  0,  no  boundary  conditions  should  be  given  there. 
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For  a  inulUclirncnsional  problem  like  this,  (  becomes  a  wave  number  vector 
and  the  group  speed  ^.7)  generalizes  to  a  vector  group  velocity  given  by 


C  = 


(3.6) 


where  denotes  the  gradient  with  respect  to  As  in  the  one-dimensional  case, 
diiTcrence  schemes  not  only  lead  to  incorrect  group  speeds,  but  may  cause  propagation 
in  the  wrong  direction — which  means  at  any  angle  in  the  (x,y)  plane  whatsoever. 
See  [13]  for  a  discussion  with  examples.  In  two  dimensions,  Thm.  I  becomes:  if  a 
difference  model  of  (S.5j  admits  a  solution  consisting  of  waves  with  group  velocity  C 
pointing  into  x  >  0  (i.e.  with  ^  0^,  is  unstable.  If  one  such  wave  has  >  0, 
then  unbounded  growth  in  will  take  place,  as  in  Thm.  2. 

For  example,  suppose  (3.5)  is  modeled  by  the  leap  frog  formula 


—  V, 


n-l 


(3.7) 


The  dispersion  relation  for  this  scheme  is 


sin  wk  =  —X  sin  ^h  —  \  sin  77/1, 


where  {  s  ((r’7)>  &nd  from  (3.6)  there  follow  the  group  velocity  components 


cos^h  ^  ^  cosv^h 
cos  ujk'  ^  coswA* 


As  usual,  these  reduce  to  the  ideal  value  C  =  (-1,-1)  for  f/i,  ujk  0.  If  we  look 
at  parasites,  on  the  other  hand,  we  see  that  a  sawtooth  form  in  x  or  y  negates  C,  or 
Cy,  respectively,  and  a  sawtooth  in  t  negates  both.  One  has 


ih,rih,  wk 

X 

(a) 

(0,0,0),  ()r,i,w) 

(-1,-1) 

(b) 

(;r,0,0),  (0,ir,ir) 

(•♦■1,-1) 

(c) 

(0,7r,0),  (x,0,7r) 

(-1,+1) 

(d) 

(ir,7r,0),  (0,0,  tt) 

(+1,+1) 

Thus  parasites  can  travel  in  any  of  the  directions  at  45°  to  the  grid.  If  any  parasite 
of  form  (b)  or  (d)  is  permitted  by  the  boundary  conditions,  the  dilTerence  model  is 
unstable. 

Abarbanel  et  al.  consider  various  boundary  formulas.  Four  of  these  arc  space 
extrapolation  and  skewed  space  extrapolation, 

S:  (£,- l)«t;5+'  =0, 

SS:  (E,Ey  -  l)*v5^'  =  0, 

and  spaee/time  extrapolation  and  skewed  space/time  extrapolation, 

ST:  (fi,£^r'  -  l)’wo'^'  =  0, 
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SST:  [E^EyE'C^  -  =  0. 

Here  ^yi  2ind  E^  denote  the  shift  operators  in  z,  y,  and  f.  By  counting  sign 
changes,  one  can  see  which  boundary  formulas  permit  which  sawtooths.  One  finds 


stable  sawtooths 

unstable  sawtooths 

s 

(0,0,0),(0,ir,0) 

(0,0,  t),(0,w,x) 

SS 

(0, 0, 0),(ir,ir,») 

(0, 0,  w),  (ir,  IT,  0) 

ST 

(0, 0, 0),  (0,  IT,  0),  (jr,  0,  tt), 

(tt  ,-»,») 

SST 

(0, 0, 0),  (t,  0,  t) 

(0,  jr,  jr),  {ir,n,0} 

Thus  S,  SS,  and  SST  are  all  unstable  with  LF.  It  turns  out  that  ST,  which  we  see  has 
no  sawtooth  instabilities,  is  indeed  stable. 

Other  difference  formulas  typically  permit  fewer  sawtooths,  hence  are  stable  with 
more  boundary  conditions.  Let  us  generalize  to  d  space  dimensions.  If  k  and  j  are 
d'vectors,  will  denote 

Defn.  Let  Q  be  a  scalar  difference  formula  in  d  space  dimensions.  Suppose  that 
whenever  Q  admits  a  solution  Vy  =  with  \z\  =  jic^l  =  1  for  each  i,  ki  =  I  for 
some  /,  and  group  velocity  0  €  then  it  also  admits  the  solution  Vj  =  (— 
and  this  wave  has  group  velocity  C'  G  satisfying  CJ  =  C»  for  i  I  and  CiCj  < 
0.  Then  Q  is  Zj-reversing.  Suppose  that  whenever  Q  admits  a  solution  Vy  =  hP 
with  \Ki\  =  1  for  all  t  and  group  velocity  C  €  then  it  also  admits  the  solution 
vj  =  1)”,  with  group  velocity  satisfying  CiC[  <  0  for  1  <  i  <  d.  Then 

Q  is  f-reversing. 

Now  let  Q  be  a  difference  model  of 

d 

j-i 

on  t,zi  >  0,  Xj  G  (—00,00)  for  2  <  j  <  d,  and  let  the  boundary  conditions  S,  SS, 
ST,  SST  be  extended  in  the  obvious  way.  By  the  same  arguments  as  above  we  obtain 
the  following  theorem: 

Theorem  7.  The  following  assertions  hold  in  the  stated  direction  only;  their 
converses  are  not  in  gefteral  valid. 

(i)  The  model  S,Q  is  unstable  if  Q  is  t-reversing. 

(ii)  The  model  SS,Q  is  unstable  if  Q  is  t-reversing  or  if  Q  is  Xi-reversing  and 
also  Xj-reversing  for  at  least  one  j  >  2. 

The  model  SST,Q  is  unstable  if  Q  is  Xi-reversing  and/or  t-reversing,  and 
also  Xj-reversing  for  at  least  one  J  >  2.  | 

Among  the  formulas  Q  considered  by  Abarbanel  et  al.  arc  multidimensional  versions 
of  LF,  CN,  BE,  and  MC  (MacCormack*s  scheme).  One  secs  readily  that  LF  is 
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reversing  and  z^-reversing  for  each  j,  CN  and  BE  are  Xj-rcversing  for  each  j  but  not 
f-reversing,  and  MC  is  not  reversing  in  any  variable.  It  turns  out  that  all  combinations 
of  these  schemes  with  S,  SS,  ST,  or  SST  that  arc  not  ruled  unstable  by  Thm.  7  are 
in  fact  stable. 
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